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medium by solving the integro-differential equation of the radiative transfer of res- 
onant photons. To solve the differential equations numerically we use the Weighted 
Essentially Non- oscillatory method (WENO). Although the effects of dust on radia- 
tive transfer is well known, the resonant scattering of Lya photons makes the problem 
non-trivial. For instance, if the medium has the optical depth of dust absorption and 
scattering to be t„ » 1, t » 1, and r » r^, the effective absorption optical depth 



p ,! in a random walk scenario would be equal to y/rjT^Tr). We show, however, that 

Q ! for a resonant scattering at frequency vo, the effective absorption optical depth would 

^ ! be even larger than t(vo). If the cross section of dust scattering and absorption is 

ci . frequency-independent, the double-peaked structure of the frequency profile given by 

the resonant scattering is basically dust-independent. That is, dust causes neither nar- 
^ , rowing nor widening of the width of the double peaked profile. One more result is that 

'nI" i the time scales of the Lyor photon transfer in the optically thick halo are also basically 

QQ ! independent of the dust scattering, even when the scattering is anisotropic. This is be- 

cause those time scales are mainly determined by the transfer in the frequency space, 
while dust scattering, either isotropic or anisotropic, does not affect the behavior of 
the transfer in the frequency space when the cross section of scattering is wavelength- 
independent. This result does not support the speculation that dust will lead to the 
smoothing of the brightness distribution of Lya photon source with optical thick halo. 
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1. Introduction 

hya photons have been widely applied to study the physics of luminous objects at various 
epochs of the universe, such as Lya emitters, Lya blob, damped Lya system, hya forest, fluo- 
rescent hya emission, star-forming galaxies, quasars at high redshifts as well as optical afterglow 
of gamma ray bursts (Haiman et al. 2000; Fardal et al. 2001; Dijkstra & Loeb. 2009; Latif et 
al. 2011). The resonant scattering of Lyo- photons with neutral hydrogen atoms has a profound 
effect on the time, space and frequency dependencies of hya photons transfer in an optically thick 
medium, hya photons emergent from an optically thick medium would carry rich information of 
photon sources and halo surrounding the source of the hya photon. The profiles of the emission 
and absorption of the hya radiation are powerful tools to constrain the mass density, velocity, tem- 
perature and the fraction of neutral hydrogen of the optically thick medium. Radiation transfer of 
hya photons in an optically thick medium is fundamentally important. 

The radiative transfer of hya photons in a medium consisting of neutral hydrogen atoms has 
been extensively studied either analytically or numerically. Yet, there have been relatively few 
results which are directly based on the solutions of the integro-differential equation of the resonant 
radiative transfer. Besides the Field solution (Field 1959, Rybicki & Dell' Antonio 1994), analyt- 
ical solutions with and without dust mostly are based on the Fokker-Planck (P-F) approximation 
(Harrington 1973, Neufeld 1990, Dijkstra et al. 2006). The P-F equation might miss the detailed 
balance relationship of resonant scattering (Rybicki 2006), and therefore, the analytical solutions 
cannot describe the formation and evolution of the Wouthuysen-Field (W-F) local thermalization 
of the hya photon frequency distribution (Wouthuysen 1952, Field 1958), which is important for 
the emission and absorption of the hydrogen 21 cm line (e.g. Fang 2009). The features of the hya 
photon transfer related to the W-F local thermalization are also missed. An early effort (Adams et 
al. 1971) trying to directly solve the integro-differential equation of the resonant radiative transfer 
with numerical method. It still is, however, of a time-independent approximation. 

Recently, a state-of-the-art numerical method has been introduced to solve the integro-differential 
equation of the radiative transfer with resonant scattering (Qiu et al. 2006, 2007, 2008, Roy et al. 
2009a). The solver is based on the weighted essentially non-oscillatory (WENO) scheme (Jiang & 
Shu 1996). With the WENO solver, many physical features of the transfer of Lya photons in an 
optically thick medium (Roy et al. 2009b, 2009c, 2010), which are missed in the Fokker-Planck 
equation approximations, have been revealed. For instance, the WENO solution shows that the 
time scale of the formation of the W-F local thermal equilibrium actually is only about a few hun- 
dred times of the resonant scattering. It also shows that the double peaked frequency profile of the 
Lya photon emergent from an optically thick medium does not follow the time-independent solu- 
tions of the P-F equation. These results directly indicate the needs of re- visiting problems which 
have been studied only via the F-P time-independent approximation. 



We will investigate, in this paper, the effects of the dust on the Lya photons transfer in an 
optically thick medium. Dust can be produced at epochs of low and moderate redshifts, and even 
at redshift as high as 6 (Stratta et al. 2007). Absorption and scattering of dust have been used 
to explain the observations on hya emission and absorption (Hummer & Kunasz 1980), such as 
the escaping fraction of Lya photons (Hayes et al. 2010, 2011, Blanc et al. 2010); the redshift- 
dependence of the ratio between Lya emitters and Lyman Break galaxies (Verhamme et al. 2008); 
and the "evolution" of the double-peaked profile (Laursen et al. 2009). 

Nevertheless, it is still unclear whether the time scale of photon escaping from optically thick 
halo will be increasing (or decreasing) when the halo is dusty. It is also unclear whether the effects 
of dust absorption can be estimated by the random walk picture (Hansen & Oh 2006). As for the 
dust effect on the double-peaked profile, the current results given by different studies seem to be 
contradictory: some claims that the dust absorption leads to the narrowing of the double-peaked 
profile (Lauresen et al 2009), while others result that the width between the two peaks apparently 
should be increasing due to the dust absorption (Verhamme et al. 2006). We will focus on these 
basic problems, and examine them with the solution of the integro-differential equation of radiative 
transfer. 

This paper is organized in the following way: section 2 presents the theory of the Lya photon 
transfer in an optically thick medium with dust. The equations of the intensity and flux of resonant 
photons in a dusty medium are given. We will study three models of the interaction between dust 
and photons: (1) dust causes only scattering with photons; (2) dust causes both scattering and 
absorption; and (3) dust causes only absorption of photons. Section 3 gives the solutions of Lyo- 
photons escaping from an optically thick spherical halos with dust. The dusty effect on the double- 
peaked profile will be studied in Section 4. The discussion and conclusion are given in Section 5. 
Some mathematical derivations of the equations and numerical implementation details are given in 
the Appendix. 



2. Basic theory 

2.1. Radiative transfer equation of dusty halo 

We study the transfer of Lyo- photons in a spherical halo with radius R around an optical 
source. The halo is assumed to consist of uniformly distributed HI gas and dust. The optical depth 
of HI scattering over a light path dl is dr = cr{v)n}iidl, where hhi is the number density of HI, and 
(t(v) is the cross section of the resonant scattering of Lya photons by neutral hydrogen, which is 



given by 

where (f)(x,a) is the normalized Voigt profile (Hummer 1965). As usual, the photon frequency y 
in eq.(l) is described by the dimensionless frequency x = (v - vq)/Avd, with vq = 2.46 x 10^^ s^^ 
being the resonant frequency, Avo = vo{vt/c) = 1.06 x 10^'(r/10"^)^''^ Hz the Doppler broadening, 
Vt = ^jlksT jm the thermal velocity, and T the gas temperature of the halo. (To/tt^^^ is the cross 
section of scattering at the the resonant frequency vq. The parameter a in eq.(l) is the ratio of 
the natural to the Doppler broadening. For the hya line, a = 4.7 x 10 '^(r/10"^)^^''^. The optical 
depth of Lya photons with respect to HI resonant scattering is Ts{x) = nmRc(x) = TQ(p{x, a), where 
To = nmo-oR. 

If the absorption and scattering of dust are described by effective cross-section per hydrogen 
atom (Tdix), the total optical depth is given by 

t{x) = To(p{x, a) + Tdix) (2) 

where the dust optical depth Td{x) = nuicrd(x)R. This is equal to assume that dust is uniformly 
distributed in IGM. The effects of inhomogeneous density distributions of dust (Neufeld 1991; 
Haiman & Spaans 1999) will not be studied in this paper. 

The radiative transfer equation of Lya photons in a spherical halo with dust is given by 

dl dl {\-i/)dI dl _ 
di] dr r dfj. dx 

-(l)(x;a)I + I 'R(x,x';a)I(r],r,x',fi')dx'diJ.'/2 (3) 

-k(x)I + Ak(x) I 'R'^ix, x';ix,ix;a)I{ri, r, x' , fi')dx' dfx' + S 

where I{t,rp,x,iu) is the specific intensity, which is a function of time t, radial coordinate r^, 
frequency x and the direction angle, fj. = cos 6, with respect to the radial vector r. 

In eq.(3), we use the dimensionless time rj defined asrj = crimo-Qt and the dimensionless radial 
coordinate r defined as r = nmcrorp. That is, r] and r are, respectively, in the units of mean free 
flight-time and mean free path of photon vq with respect to the resonant scattering without dust 
scattering and absorption. Without resonant scattering, a signal propagates in the radial direction 
with the speed of light, the orbit of the signal is then r = r] + const. With dimensionless variable, 
the size of the halo R is equal to tq. 

The re-distribution function 'R(x, x';a) gives the probability of a photon absorbed at the fre- 
quency x', and re-emitted at the frequency x. It depends on the details of the scattering (Henyey & 



Greestein 1941; Hummer 1962; Hummer 1969). If we consider coherent scattering without recoil, 
the re-distribution function with the Voigt profile can be written as, 

n{x,x'\a)= (4) 
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where x^^-a = min(;c, x') and x^ax = niax{x, x'). In the case of a = 0, i.e. considering only the 
Doppler broadening, the re-distribution function is 

«(x,.v') = ierfc[max(W,|x-|)l. (5) 

Xco I/O 2 

^Tiix, x')dx' = (p(x,0) = n~ ' e~^ . 

With this normalization, the total number of photons is conserved in the evolution described by 

equation (3). That is, the destruction processes of hya photons, such as the two-photon process 

(Spitzer & Greenstein 1951; Osterbrock 1962), are ignored in equation (3). The recoil of atoms is 

also not considered in equation (4) or (5). The effect of recoil actually is under control (Roy et al. 

2009c, 2010). We will address it in next section. 

The absorption and scattering of dust are described by the term k{x)I of eq.(3), where k{x) = 
crjo-o, which is of the order of 10^^(7/ lO^)!/^ (pj-aine & Lee 1984; Draine 2003). The term with 
A of eq.(3) describes albedo, i.e. A = crjo-d, where cr^ is the cross section of dust scattering. 
Generally, A lies approximately between 0.3 and 0.4 (Pei 1992; Weingartner & Draine 2001). 

Since dust generally is much heavier than a single atoms, the recoil of dust particles can be 
neglected when colliding with a photon. Under this "heavy dust" approximation, photons do not 
change their frequency during the collision with dust. The redistribution function of dust 'R'' is 
independent of x and x', and is simply given by a phase function as 

'R'(M,fi') = ^ r diP' . I" {-,3/2 = Z ^^^s'Pi(^)Pi^')^ (6) 

471- Jo (1 +r-2^y")^/2 j_u 2 

where jx = ixjx' + y(l-ju^)(l-/?^cos0' and Pi is the Legendre function. The factor g in eq.(6) is 
the asymmetry parameter. For isotropic scattering, ^ = 0. The cases of ^ = -i-l and -1 correspond 
to complete forward and backward scattering, respectively. Generally, the factor ^ is a function 
of the wavelength. For the Lya photon, we will take g = 0.73 for realistic dust scattering (Li & 
Draine 2001). The integral of eq.(6) is performed in Appendix A. 

In eq. (3), the term with the parameter y is due to the expansion of the universe. If n^ is 
equal to the mean of the number density of cosmic hydrogen, we have y = r^^p, and tgp is the 
Gunn-Peterson optical depth. Since the Gunn-Peterson optical depth is of the order of 10^ at high 



redshift (e.g. Roy et al. 2009c), the parameter y is of the order of 10 ^ - 10 ^. Therefore, if the 
optical depth of halos is equal to or less than 10^, the term with y of eq.(3) can be ignored. 

In eq.(3) we neglect the effect of collision transition from H(2p) state to H(2s) state, which 
can significantly affect on the escape of hya photons when HI column density is higher than 10^^ 
cm"^ and dust absorption is very small (Neufeld, 1990). This generally is out of the parameter 
range used below. We are also not considering the effects of bulk motion of the medium of halos 
(e.g. Spaans & Silk 2006, Xu & Wu, 2010). 



2.2. Eddington approximation 

Eq.(6) indicates that the transfer equation (3) can be solved with the Legendre expansion 
/(?7, r, x,iJ.) = YjI h{f], r, x)Pi{p). If we take only the first two terms, 1 = and 1, it is the Eddington 
approximation as 

/(?7, r, X, n) - J(7], r, x) + 3yuF(?7, r, x) (7) 

where 

1 r' 1 r' 

J(r],r,x) = - I I{r],r,x,n)dix, F(r],r,x) = - I ixl{r],r,x,ix)dn. (8) 

They are, respectively, the angularly averaged specific intensity and flux. Defining j = r^J and 
/ = r^F, Eq.(3) yields the equations of j and / as 

T- + % = -('^-A)KJ-cf>(x;a)j+ f 'R(x,x';a)jdx' + y^ + r^S, (9) 

or] or J ox 

^ + \^-¥ = -ii-Ag)Kf + y^-cf>ix;a)f. (10) 

Orj 5 or 5 r ox 

The mean intensity j(ri, r, x) describes the x photons trapped in the position r at time 77 by the 
resonant scattering, while the flux /(?7, r, x) describes the photons in transit. 

The source term S in the equations (3) and (9) can be described by a boundary condition of j 
and / at r = r^. We can take r^ = 0. Thus, the boundary condition is 

M 0, X) = 0, fin, o,x) = s ,4 six), (11) 

where 5o, and (psix) are, respectively, the intensity and normalized frequency profile of the sources. 
Since equation (3) is linear, the solutions of jix) and fix) for given 5o = 5 are equal to S jiix) 
and S fiix), where jiix) and fiix) are the solutions of 5o = 1- On the other hand, the equation (3) 
is not linear with respect to the function (psix). The solution fix) for a given (psix) is not equal to 
(f>six)fiix), where fiix) is the solution of (psix) = 1. 



In the range outside the halo, r > R, no photons propagate in the direction ju < 0. The 
boundary condition at r = R given by J fj.J(r],R, x,iu)dfj. = is then (Unno 1955) 

MR,x) = 2f(u,R,x). (12) 

There is no photon in the field before t = 0. Therefore, the initial condition is 

j(0,r,x) = f{0,r,x) = 0. (13) 

We will solve equations (9) and (10) with boundary and initial conditions eqs.(ll) - (13) by using 
the WENO solver (Roy et al. 2009a, b, c, 2010). Some details of this method is given in Appendix 
B. 



2.3. Dust models 

We consider three models of the dust as follows: 

I. pure scattering, A = 1, ^ = 0.73: dust causes only anisotropic scattering, but no absorption; 

II. scattering and absorption. A = 0.32, g = 0.73: dust causes both absorption and anisotropic 

scattering. 

III. pure absorption. A = 0: dust causes only absorption, but no scattering; 

Models I and III do not occur in reality. They are, however, helpful to reveal the effects of pure 
scattering and absorption on the radiative transfer. 

Since k{x) is on the order of 10"^, its effect will be significant only for halos with halos with 
optical depth tq > 10^, and ignorable for tq < 10^. To illustrate the dust effect, we use halos 
of R = To < 10"*, and take larger k to he ^ 10 "* - 10"^. We also assume that k is frequency- 
independent. We consider below only the case of grey dust, i.e. k is independent of frequency x. 
This certainly is not realistic dust. Yet, the frequency range given in solution below mostly are in 
the range \x\ < 4. Therefore, the approximation of grey dust would be proper if cross section of 
dust is not significantly frequency dependent in the range \x\ < 4. 



2.4. Numerical example: Wouthuysen-Field thermalization 

As the first example of numerical solutions, we show the Wouthuysen-Field (W-F) effect, 
which requires that the distribution of Lyor photons in the frequency space should be thermalized 



near the resonant frequency vq. The W-F effect illustrates the difference between the analytical 
solutions of the Fokker- Planck approximation and that of eqs. (9) and (10). The former can not 
show the local thermalization (Neufeld 1990), while the latter can (Roy et al. 2009b). All problems 
related to the W-F local thermal equilibrium should be studied with the integro-differential equation 

(3). 
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Fig. 1. — The mean intensity j(r], r, x) at rj = 500 and r = 100 for dust models I (left panel), II 
(middle panel) and III (right panel). The source is So = 1 and cpX^) = (1/ ^ln)e'^ . The parameter 
a = 10"^. In each panel, k is taken to be 0, 10""^, 10"^ and 10"^. 

Figure 1 presents a solution of mean intensity jirj, r, x) at time radial rj = 500 coordinate 
r = 10^ for halo with size i? » r = 10^. The three panels correspond to dust models I (left 
panel), II (middle panel) and III (right panel). The source is taken to have a Gaussian profile 
^s(x) = (1/ ^ln)e~'^ and unit intensity 5o = 1. The solutions of Figure 1 actually are independent 
of R, li R » 10^. The intensity of j is decreasing from left to right in Figure 1, because the 
absorption is increasing with the models from I to III. 

A remarkable feature shown in Figure 1 is that all j{ri, r, x) have a flat plateau in the range 
\x\ < 2. This gives the frequency range of the W-F local thermalization (Roy et al, 2009b, c). The 
range of the fiat plateau \x\ < 2 is almost dust-independent, either for model I or for models II and 
III. This is expected, as neither the absorption nor scattering given by the k term of eq.(3) changes 
the frequency distribution of photons. The redistribution function (6) also does not change the 
frequency distribution of photons. This point can also be seen from eqs. (9) and (10), in which the 
K terms are frequency-independent. The evolution of the frequency distribution of photons is due 
only to the resonant scattering. 

Since thermalization will erase all frequency features within the range \x\ < 2, the double- 
peaked structure does not retain information of the photon frequency distribution within \x\ < 2 at 
the source. That is, the results in Figure 1 will hold for any source So((>s{x) with arbitrary cpsix) 
which is non-zero within \x\ < 2 (Roy et al. 2009b, c). This property can also be used as a test of 
the simulation code. It requires that simulation results of flat plateau should be hold, regardless of 



the source to be monochromatic or with finite width around vq. 



3. Dust effects on photon escape 



3.1. Model I: scattering of dust 

To study the effects of dust scattering on the Lya photon escape, we show in Figure 2 the flux 
/(?7, r, x) of Lyor photons emergent from halos at the boundary r = R = 10^ for Model I. The three 
panels of Figure 2 correspond to a: = 10 "^, 10"^, and 10"^ from left to right, respectively. The 
source starts to emit photons at ?/ = with a stable luminosity 5o = 1, and with a Gaussian profile 

0,W = (1/V^K^. 




Fig. 2. — Flux f(7], r, x) of Lya photons emergent from halos at the boundary R = 10^, and for the 
dust model I A = 1, ^ = 0.73. The parameter k is taken to be 10^"^ (left), 10"^ (middle) and 10"^ 
(right). The source is So = 1 and ^s(x) = (1/ ^fn)e''' . The parameter a = lO"-'. 



Figure 2 clearly shows that the time-evolution of /(?/, r, x) is /c-independent. Although the 



cross section of dust scattering increases about 100 times from k 
the left and right panels in Figure 2 actually are almost identical. 



10 to /(• = 10 , the curves of 



According to the scenario of "single longest excursion", photon escape is not a process of 
Brownian random walk in the spatial space, but a transfer in the frequency space (Osterbrock 
1962; Avery & House 1968; Adams, 1972, 1975; Harrington 1973; Bonilha et al. 1979). Photon 
will escape, once its frequency is transferred from \x\ < 2 to \x\ > 2, on which the medium is 
transparent. On the other hand, dust scattering given by the redistribution function eq.(6) does not 
change photon frequency. Dust scattering has no effect on the transfer in the frequency space. 

Moreover, photons with frequency \x\ < 2 are quickly thermalized after a few hundred res- 
onant scattering. In the local thermal equilibrium state, the angular distribution of photons is 
isotropic. Thus, even if the dust scattering is anisotropic ^ ^ with respect to the direction of the 



10 



incident particle, the angular distribution will keep isotropic undergoing ag i^ scattering. Hence, 
dust scattering also has no effect on the angular distribution. 



3.2. Model III: absorption of dust 



Similar to Figure 2, we present in Figure 3 the flux of Model III, i.e. dust causes only ab- 
sorption without scattering. All other parameters of Figure 3 are the same as in Figure 2. In the 
left panel of Figure 3, the curves at the time i] = 2000 and 3000 are the same. It means the flux 
/(?7, R, x) at the boundary R is already stable, or saturated at the time r] > 2000. The small dif- 
ference between the curves of rj = 1000 and r] > 2000 of the left panel indicates that the flux is 
still not yet completely saturated at the time r] = 1000. However, comparing the middle and right 
panels of Figure 3, we see that for k = 10"^, the flux has already saturated at 77 = 1600, while it 
has saturated at 77 = 800 for k = 10^. That is, the stronger the dust absorption, the shorter the 
saturation time scale. The time scales of escape or saturation do not increase by dust absorption, 
and even decrease with respect to the medium without dust. Stronger absorption leads to shorter 
time scale of saturation. 
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Fig. 3. — Flux f{ri,r,x) of hya photons emergent from halos at the boundary r = R = 10^. 
The parameters of the dust are A = and k = 10^"^ (left), 10"^ (middle) and 10^^ (right). Other 
parameters are the same as in Figure 2. 



Obviously, dust absorption does not help in producing photons for the "single longest excur- 
sion". Therefore, dust absorption can not make the time scale of producing photons for "single 
longest excursion" to be smaller. However, dust absorptions are efli'ective in reducing the number 
of photons trapped in the state of local thermalized equilibrium \x\ < 2 (see also §4.2). This leads 
to the fact that the higher the value of k, shorter the time scale of saturation. 
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3.3. Effective absorption optical depth 

Since hya photons underwent a large number of resonant scattering before escaping from 
halo with optical depth tq » 1, it is generally believed that a small absorption of dust will lead to a 
significant decrease of the flux. However, it is still unclear what the exact relationship between the 
dust absorption and the resonant scattering is. This problem should be measured by the eff'ective 
optical depth of dust absorption of Lya photons in 7? = tq » 1 halos. 

To calculate the efli'ective optical depth, we first give the total flux of Lya photons emergent 
from halo of radius R, which is defined as F(r]) = J /(?], R, x)dx. Figure 4 plots F{if) as a function 
of time ?7 for halo with sizes R = tq = \(fi and 10"^. The curves typically are the time-evolution of 
growing and then saturating. The three panels correspond to the dust models I, II and III from left 
to right. The upper panels are of 7? = 10^, and lower panels for i? = 10"*. In each panel of 7? = 10^, 
we have three curves corresponding to /c = 10""^, lO""* and 10"^, respectively. In cases oiR = \0^, 
we take k = 10""^ and lO^-'. 
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Fig. 4. — The time evolution of the total flux F{ri) at the boundary of halos with R = tq = 10^ 
(upper panels), and 7? = tq = lO'* (lower Panels). The source of So = 1 and (ps(x) = (1/ V^)e 
starts to emit photons at time rj = 0. The parameters of dust are (A = 1,^ = 0.73) (left); (A = 
0.32, g = 0.73) (middle) and A = (right). In each panel of R = 10^ k is taken to be 10 \ 10^^ 
and 10^2. In the cases ofR= lOl k is taken to be 10"\ 10 I 



12 



Theleftpanelof Figure 4 shows that the three curves of a: = 10^, 10^ and 10^^ are almost the 
same. This is consistent with Figure 2 that for Model I, the time-evolution of / are /c-independent 
for the pure scattering dust. For the pure absorption dust (the right panel of Figure 4), the saturated 
flux is smaller for larger k. We can also see from Figure 4 that the time scale of approaching 
saturation is smaller for larger k. The result of model II is in between that for models I and III. 

With the saturated flux of Figure 4, one can define the effective absorption optical depth by 
^effect = -(^/k) ln/5. The results are shown in Table 1, in which r^ is the dust absorption depth. It 
is interested to see that the eff"ective absorption optical depth is always equal to about a few times 
of the optical depth of resonant scattering tq, regardless whether r^ is less than 1. Namely, the 
effective absorption depth reflect of dust is roughly proportional to tq. 

Table 1 . Effective absorption optical depth reflect 





Model II 


Model III 


R = To 


K 


Ta 


fs 


Teflect 


Ta 


fs 


Teffect 


102 
102 
102 
104 
104 


10-4 
10-3 
10-2 
10-4 
10-3 


0.0068 

0.068 

0.68 

0.68 

6.8 


0.978 

0.760 

0.116 

6.28 X 10-2 

4.07 X 10-^ 


2.2 X 102 
2.7 X 102 
2.2 X 102 
2.8x104 
1.5x104 


0.01 
0.10 
1.00 
1.00 
10.0 


0.963 

0.670 

0.057 

3.02 X 10-2 

2.87 X 10-'^ 


3.8x102 
4.0 X 102 
2.9 X 102 
3.5 X 104 
1.97x104 



According to the random walk scenario, if a medium has optical depths of absorption r^ and 
scattering r,, the effective absorption optical depth should be equal to reflect = V'^aC'^a + ^s) (Ry- 
bicki & Lightman 1979). However, the results of the last line of Table 1 show that the random walk 
scenario does not work for the dust effect on resonant photon transfer. This result is consistent with 
Figures 2 and 3. When optical depth of dust is lower than the optical depth of resonant scattering 
To, the time scale of photon escaping basically is not affected by the dust, but is proportional to Tq, 
and therefore, the absorption is also proportional to tq. 



3.4. Escape coefficient 



With the total flux, we can define the escaping coefficient of Lya photon as fesdn^ "^o) = 
F{ri)/FQ, where Fq is the flux of the center source. Figure 5 shows fesdri. To) at three times 7/ = 
5 X 103, 104 and 3.2xl04 for Model II and k = 10^. At 77 = 5 x 103, the flux of halos with tq < 103 
is saturated. At 77 = 104, halos with ro < 3 x 103 are saturated, and all halos of tq < 104 are 
saturated at 77 = 3.2 x 104. 
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Fig. 5. — Escaping coefficient /esc(?7) as a function of the optical depth To of halo at time ?7 = 5x10^, 
\Q\ and 3.2x10^ from bottom to up. Dust is modeled by II, A = 0.32, g = 0.73, and k = 10~l 

4. Dust effects on double-peaked profile 



4.1. Dust and the frequency of double peaks 

A remarkable feature of Lya photon emergent from optically thick medium is the double- 
peaked profile. Figures 1, 2 and 3 have shown that the double peak frequencies x+ = |x_| are 
almost independent of either the scattering or the absorption of dust. In this section, we consider 
halos with size R or tq larger than 10^. Figure 6 presents the double peak frequency \x+\ as a 
function of utq, where the parameter a is taken to be 10"^ (left) and 5 x 10"^ (right). Comparing 
the curves with dust and without dust in Figure 6 we can say that the dust effect on \x+\ is very 
small till flT() = aR = 10^. 

In the range aTo < 20, the |jc±|-to relation is almost flat with \x+\ - 2. It is because the double- 
peaked profile is given by the frequency range of the locally thermal equilibrium. The positions 
of the two peaks, x+ and x^, basically are at the maximum and minimum frequencies of the local 
thermalization. The frequency range of the local thermal equilibrium state is mainly determined by 
the Doppler broadening, and weakly dependent on To. Thus, we always have x± - ±2. When the 
optical depth is larger, utq ~ 10^, more and more photons of the flux are attributed to the resonant 
scattering by the Lorentzian wing of the Voigt profile. In this phase, \x+\ will increase with tq. 

Figure 6 shows also a line x+ = +(aro)^^^, which is given by the analytical solution of the 
Fokker-Planck approximation, in which the Doppler broadening core in the Voigt profile is ignored 
(Harrington 1973, Neufeld 1990, Dijkstra 2006). The numerical solutions of eqs (3) or (9) and 
(10) deviate from the (aToY^^-law at all parameter range of Figure 6. The deviation at utq < 20 
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Fig. 6. — The two-peak frequencies x^ = |x_| as a function of aro. The parameter a is taken to be 
10"^ (left) and 5 x 10"^ (right). Dust model III (pure absorption) is used, and k is taken to be 10"^. 
The dashed straight line gives log .x:+-log ar with slope 1/3, which is to show the (aT)^^^-law of x±. 

is due to that the Doppler broadening core in the Voigt profile is ignored in the Fokker-Planck 
approximation, and then, no locally thermal equilibrium can be reached. Therefore, in the range 
aTo < 20, \x+\ of the WENO solution is larger than the (aToY^^-law. In the range of utq > 20, 
the Fokker-Planck approximation yields a faster diffusion of photons in the frequency space. This 
point can be seen in the comparison between a Fokker-Planck solution with Field's analytical 
solution (Figure 1 in Rybicki & Dell' Antonio 1994). In this range, the numerical results of \x+\ is 
less than the ((2To)^''^-law. 



4.2. No narrowing and no widening 

The dust effect has been used to explain the narrowing of the width between the two peaks 
(Laursen et al. 2009). Oppositely, it is also used to explain the widening of the width between the 
two peaks (Verhamme et al. 2006). However, Figures 1, 2, 3 and 6 already show that the width 
between the two peaks of the profile is very weakly dependent on dust scattering and absorption. 
This result supports, at least in the parameter range considered in Figures 1, 2, 3, neither the 
narrowing nor the widening of the two peaks. 

If dust absorption can cause narrowing, the absorption should be weaker at \x\ ~ 0, and 
stronger at \x\ > 2. Similarly, if dust absorption can cause widening, the absorption should be 
weaker at \x\ ~ 2, and stronger at \x\ ~ 0. To test these assumptions. Figure 7 plots ln[/(?7, r,x,K = 
0)/f{rj, r, X, k)] as a function of x. It measures the jc(frequency)-dependence of the flux ratio with 
and without dust absorption. We take large -q, and then the fluxes in Figure 7 are saturated. Figure 
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Fig. 7. — ln[/(7/, r, .)c, a: = 0)//(?7, r, .jc, /<•)] as function of x for model II (up), and III (bottom), and 
K = 10"^ (left) and 10~^ (right). Other parameters are the same as in Figure 2. 

7 shows that the absorption in the range \x\ < 2 is much stronger than that of \x\ > 2, and therefore, 
the assumption of the narrowing is ruled out. Figure 7 shows also that the curves of \n[f{ri, r,x,K = 
0)1 f{ri, r,x,K= 10"^)] are almost flat in the range \x\ < 2. Therefore, the assumption of widening 
of the two peaks can also be ruled out. 

Since the cross sections of dust absorption and scattering are assumed to be frequency- 
independent. Eqs. (9) and (10) do not contain any frequency scales other than that from resonant 
scattering. However, either narrowing or widening would require to have frequency scales different 
from that of resonant scattering. This is occurence is not possible if the dust is gray. 



4.3. Profile of absorption spectrum 

If the radiation from the sources has a continuum spectrum, the eff'ect of neutral hydrogen 
halos is to produce an absorption line at y = vq. The profile of the absorption line can also be 
found by solving equations (9) and (10), but replacing the boundary equation (1 1) by 



j{r], 0, x) = 0, 



f{r],0,x) = So. 



(14) 
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That is, we assume that the original spectrum is flat in the frequency space. The spectrum of the 
flux emergent from halo ofR= 10^ and 10"^ with central source of eq.(14) for dust models I, II and 
III are shown in Figure 8. All curves are for large rj, i.e. they are saturated. 
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Fig. 8. — The spectrum of the flux emergent from halo ofR = 10^ (upper panels) and 10"^ (lower 
panels) with central source of eq.(14) for the dust model I (left), II (middle) and III (right). Other 
parameters are the same as in Figure 2. 



The optical depths at the frequency \x\ > 4 are small, and therefore, the Eddington approx- 
imation might no longer be proper. However, those photons do not strongly involve the resonant 
scattering, and hence they do not significantly affect the solution around x = 0. The solutions of 
Figure 8 is still useful to study the profiles of/ around x = 0. 

The flux profile of Figure 8 typically are absorption lines with width given by the double 
peaks similar to the double peaked structure of the emission line. The flux at the double peaks 
is even higher than the flat wing. It is because more photons are stored in the frequency range 
\x\ < 2. According to the redistribution function eq.(4), the probability of transferring a x' photon 
to a \x\ < \x'\ photon is larger than that from \x'\ to \x\ > \x'\. Therefore, if the original spectrum 
is flat, the net effect of resonant scattering is to bring photons with frequency \x\ > 2 to \x\ < 2. 
Photons stored \x\ < 2 are thermalized, and therefore, in the range \x\ < 2, the profile will be the 
same as the emission line, and the double peaks can be higher than the wing. It makes the shoulder 
at \x\ ~ 2. 
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As expected, for model I (left panels of Figure 8), the double profile is completely A:-independent. 
Dusty scattering does not change the flux and its profile. For models II and III, the higher the k, the 
lower the flux of the wing, because the dust absorption is assumed to be frequency-independent. 
The positions of the double peaks, x, in the absorption spectrum are also A:-independent. This 
once again shows that dust absorption and scattering causes neither narrowing nor widening of the 
double-peaked profile. However, for higher k the flux of the peaks is lower. When the absorption is 
very strong, the double-peaked structure might disappear, but will never be narrowed or widened. 



5. Discussions and conclusions 

The study of dust eff'ects on radiative transfer has had a long history related to extinction. 
However, dust efli'ects on radiative transfer of resonant photons actually have not been carefully 
investigated. Existing works are mostly based on the solutions of the Fokker-Planck approxima- 
tion, or Monte Carlo simulation. These results are important. We revisited these problems with the 
WENO solver of the integro-diff"erential equation of the resonant radiative transfer, and have found 
some features which have not been addressed in previous works. These features are summarized 
as follows. 

First, the random walk picture in the physical space will no longer be available for estimating 
the eff'ective optical depth of dust absorption. For a medium with the optical depth of absorption 
and resonant scattering to be r^ » 1, t(V()) » 1 and TsCvq) » r^, the effective absorption optical 
depth is found to be almost independent of r,,, and to be equal to about a few times of Tj(vo). 

Second, dust absorption will, of course, yield the decrease of the flux of Lya photons emer- 
gent from optical thick medium. However, if the absorption cross-section of dust is frequency 
independent, the double-peaked structure of the frequency profile is basically dust-independent. 
The double-peaked structure does not get narrowed or widened by the absorption and scattering of 
dust. 

Third, the time scales of Lyo- photon transfer basically are independent of dust scattering and 
absorption. It is because those time scales are mainly determined by the kinetics in the frequency 
space, while dust does not affect the behavior of the transfer in the frequency space if the cross 
section of the dust is wavelength-independent. The local thermal equilibrium makes the anisotropic 
scattering to be ineffective on the angular distribution of photons. Dust absorption and scattering 
do not lead to the increase or decrease of the time of storing Lyo- photons in the halos. 

The differences between the time-independent solutions of the Fokker-Planck approximation, 
or Monte Carlo simulation and the WENO solution of eq.(3) is mainly related to the W-F effect. 
Therefore, all above-mentioned features can already be clearly seen with halos of tq ~ 10^, in 
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which the W-F local thermal equilibrium has been well established. 

In this context, most calculation in this paper is on holes with tq < 10^. This range of tq 
certainly is unable to describe halos with column number density of HI larger than 10'^ cm"^ (e.g. 
Roy et al. 2010). Nevertheless, the result of tq < 10^ would already be useful for studying the 
21 cm region around high-redshift sources, of which the optical depth typically is (Liu et al 2007; 
Roy, et al. 2009c). 

where /hi is the fraction of HI. All other parameters in eq. (15) is taken from the concordance 
ACDM mode. For these objects the relation between dimensionless rj and physical time t is given 

The 21 cm emission rely on the W-F effect. On the other hand, the time-scale of the evolution of 
the 21 region is short. The effect of dust on the time-scales of Lya evolution should be considered. 

We have not considered the hya photons produced by the recombination in the ionized halo. If 
the halo is optical thick, photons from the recombination will also be thermalized. The information 
of where the photon comes from will be forgotten during the thermalization. Therefore, photons 
from recombination should not show any difference from those emitted from central sources. Only 
the photons formed at the place very close to the boundary of the halo will not be thermalized, and 
may yield different behavior. 



This research is partially supported by ARO grants W911NF-08- 1-0520 and W911NF-11-1- 
0091. 



A. Integral of the phase function [eq.(6)] 

Eq.(6) can be rewritten as 

K'iM,^') = — dcp' ^ (Al) 

4^ Jo \\-g\'\i 

where I and I' are unit vector on the direction of polar angle 6 and 9' , and azimuth angle and (p' , 
respectively. That is I • I = I' • I' = 1 and I • I' = cos y = cos 6 cos 6' + sin 6 sin 6' cos(0 - 0'), and 
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yu = COS 9, //' = COS 6. We have 

d 1 



1- 



1 



and therefore, 



dg\l-glT 2^|I-^I'p/2 2^11 -^If 



l-g^ d 1 1 



/2' 



(A2) 



The expansion with Legendre functions P/(cos 7) gives 

-. CO 

-— — = ^^'PKcosr), 

and then 

_~t!|3/2 = E 2^^'^KCOS r) + E ^'P;(C0S j). 



(A3) 



II -^r 



(A4) 



(A5) 



/=o 



Since cosy = cos cos 0' + sin0sin0'cos(0 - cp'), we have the following identity for the 
Legendre function Pi(cos y) as 



m=l 



m=\ 



{l-m)\ 



Pi(cos y) = Piicos e)Pi(cos 0') + 2 V -PTicos e)P?(cos ff) cos[m(0 - 0')]- (A6) 

^—i (I + m ) 



The integral of cp' in eq.(Al) kills the second term of eq.(A6), we have then 



'>d(,, ,,' 



1 



K'{^l,^l') = —2n 



An 

1 

2 



E llg' Piicos e)Pi(cos 6') + Y^ g'Piicos e)Pi(cos 6') 



i=\ 



1=0 



(A7) 



J] 2lg'Pi(M)Pi(M') + J] g'Pi(M)Pi(M') 



1=1 



Using the orthogonal relation J Pi(ii)Pi'(pi)dfi = ^^5ir, we have 

Ro{g) = ^fdiuf dij'R'(ju,fx') = 1, 
for which only the term / = in eq.(A7) has contribution. Similarly, 

R,(g) = ^ f dfi f diu'iuR'(n,fi') = ^ f dfi f dix'ix'R%,ix') = 0, 

^2(^) = \J "^I^J diu'iuij'R'(M,lu') = |. 
These results are used in deriving eqs.(9) and (10). 



(A8) 

(A9) 
(AlO) 
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B. Numerical algorithm 

To solve Equations (9) and (10) as a system, our computational domain is (r, x) 6 [0, rmax] x 
[^left'-^'^right]' where rmax, ■'Cieft ^^^ •"•risht ^^^ chosen such that the solution vanishes to zero 
outside the boundaries. We choose mesh sizes with grid refinement tests to ensure proper numerical 
resolution. In the following, we describe numerical techniques involved in our algorithm, including 
approximations to spatial derivatives, numerical boundary condition, and time evolution. 



B.l. The WENO Algorithm: Approximations to the Spacial Derivatives 

The spacial derivative terms in Equation (9) and (10) are approximated by a fifth-order finite 
difference WENO scheme. 

We first give the WENO reconstruction procedure in approximating ^, 

dj(rf,ri,Xj) 1 . . 
d-x " A-x^^J^'^' - ^^-"'^ 

with fixed rj = if and r = r,. The numerical flux ^y+1/2 is obtained by the fifth-order WENO 
approximation in an upwind fashion, because the wind direction is fixed. Denote 

hj = j{rf,ri,Xj), j = -2,-l,---,N + 3 

with fixed n and /. The numerical flux from the WENO procedure is obtained by 

hj+i/2 = ojihf^^i^ + a)2h%i2 + (^i,h%i2 

where ^ .^in ^^ ^^ three third-order fluxes on three different stencils given by 

,a) 1. 5. 1 






:.(3) 



11. 7. 1 



Vi/2 = y^i+i - i^m + 3^7+3- 

and the nonlinear weights oj^ are given by 

OJm . Jl 



0)m - Z^ — , <^l 



where 6 is a parameter to avoid the denominator to become zero and is taken as 6 = 10^^. The 
linear weights ji are given by 

_3_ _ 3 _ J_ 



21 



and the smoothness indicators I5i are given by 



13 ,1 , 



To approximate the r-derivatives in the system of Equations (9) and (10), we need to perform 
the WENO procedure based on a characteristic decomposition. We write the left-hand side of 
Equations (9) and (10) as 

Uf + Au^, 

where u = (7, fY and 

is a constant matrix. To perform the characteristic decomposition, we first compute the eigenvalues, 
the right eigenvectors and the left eigenvectors of A and denote them by A, R and R'^. We then 
project u to the local characteristic fields v with v = R^u. Now u, + Auy of the original system is 
decoupled as two independent equations as v, + Av^. We approximate the derivative v^^ component 
by component, each with the correct upwind direction, with the WENO reconstruction procedure 
similar to the procedure described above for ^. In the end, we transform v^ back to the physical 
space by u^ = R\r. We refer the readers to Cockbum et al. 1998 for more implementation details. 



B.2. Numerical Boundary Condition 

To implement the boundary condition (12), we also need to perform a characteristic decompo- 
sition as discussed above. Using the same notation as before, we project u to the local characteristic 
fields V with v = R'^u. Denote v = (vi, vi)^, now u^ -1- Aur of the original system is decoupled to 
two independent scalar operators given by 

dvi dvi dv2 dv2 

dt dr ' dt dr 

where A\ = ^ and A2 = ~ 3 • The characteristic line starting from the boundary r = rmax for the 
first equation is pointing outside the computational domain while the one for the second equation 
is pointing inside. For well-posedness of our system, we need to impose the boundary condition 
there as 

V2 = av\ + /3 
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with constants a and /?. We can calculate the values of a and yS based on equation (12) and the left 
and right eigenvectors of A. For example, if we take 



R 



( 2H 2H \ 

2 2 

i _i 

V 2 2 J 



we can compute that a = 7 - 4y/3 and /3 = 0. We use extrapolation to obtain the value of vi and 
then compute the value V2. In the end, we transfer v back to the physical space by u = R\. 



B.3. Time Evolution 

To evolve in time, we use the third-order TVD Runge-Kutta time discretization (Shu & Osher 
1988). For system of ODEs m, = L(u), the third order Runge-Kutta method is 

u^^'> = u" + ATLiu",T"), 

m(2) = L" + ^iu^'^ + ATLiu^'\T" + At)), 

m"+i = iu" + \{vP + AtL(u^^\t" + I At)). 
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